function [handles] = ProcessImageGUI(handles)

% [handles] = ProcessImageGUI(handles)
% Top-level image analysis function for our experiment.
% Input:
    % handles - A structure that must contain at least the following values:
        % handles.do2DFit - boolean
        % handles.dpath - string, full path to file
        % handles.autoROI - boolean
        % handles.do1DFit - boolean
        % handles.ROI - [1 x 4] matrix of integers
        % handles.useBackground - boolean
        % handles.scaleBackground - boolean
        % If it contains handles.figure1, it is assumed the function is
        % called from the GUI and will act accordingly (try to update GUI
        % display).
% Returns:
    % handles - Same as the input structure but with the following
    % values added:
        % handles.numAtoms
        % handles.pX
        % handles.pY
% Note 'handles' can just be a path to the HDF file, in which case the
% function assumes default values for the handles structure. 

% TODO:
% handles.total_num_counts
%   Save to file:
%       Integrals of x- and y- data (but over which ROI?)
%       num of atoms
%       fit results
%   Choose fitting function
%   if 'doFitData' is 0, then don't display flat lines

if isa(handles,'char')
    hdfFilePath = handles;
    
    fit2d = 0;
    autoROI = 0;
    doFitData = 1;
    userROI = [1 1392 1 1040];
    useBackground = 0;
    scaleBackground = 0;
    doDisplayImages = 1;
    saveOD = 0;
    
    % delete handles so it can be over-written later as a struct with no
    % warning from MATLAB.
%     clear handles
else
    fit2d = handles.do2DFit;
    hdfFilePath = handles.dpath;
    autoROI = handles.autoROI;
    doFitData = handles.do1DFit;
    userROI = handles.ROI;
    useBackground = handles.useBackground;
    scaleBackground = handles.scaleBackground;
    doDisplayImages = 0;
    saveOD = handles.saveOD;
end

% To check if function has been called from GUI or not,
% if it has, update the GUI appropriately, if not, just
% analyse the data and pass the results out in the structure.
% TODO: save analysis results in HDF file?
calledFromGUI = isfield(handles,'figure1');

% Inputs from hdf5 file
rawimages = cast(h5read(hdfFilePath,'/Image/image'),'double')/4; % Divide by 4 as 14-bit data from camera is multiplied by 4 (presumably when cast as a 16-bit integer).
img_var = h5read(hdfFilePath,'/Image/configuration');
img_time = h5read(hdfFilePath,'/Image/exposuretime');
type_img = h5read(hdfFilePath,'/Image/typeimaging');
CCDbin = 1; % Amount of binning at camera sensor; '1' is 1x1 binning and '2' is 2x2 binning

%Software binning for 2D fits
bindim=[4,4];  % Software binning size

% Libraries used in this function
expp = ParamImgSys(img_var);
cons = Constants();

% Size of hdf5 image group
[rows,cols,files] = size(rawimages); % rows = horiz. pixels; cols = vert. pixels

%%% TODO: Check validity of image, has pixels, correct size, has been
%%% binned etc.
if rows*cols == 0
    error('No image in HDF file (returned empty matrix)')
end

ovflw = 'no';   % setting 'ovflw' flag to 'no'
if  type_img == 3 % Absorption
    if useBackground
        elec = rawimages(:,:,3);
    else
        elec = 0;
    end
    
    if scaleBackground
        sampleIndX = 1:300; % Atoms fall in +x direction so there should never be atoms in this region. 
        sampleIndY = 1:300;
        scaleS = mean(mean(rawimages(sampleIndX,sampleIndY,1)));
        scaleB = mean(mean(rawimages(sampleIndX,sampleIndY,2)));
    else
        scaleS = 1;
        scaleB = 1;
    end
    
    num = rawimages(:,:,1) - elec;
    den = rawimages(:,:,2)*scaleS/scaleB - elec;
    
    maxcntN = max(max(num));
    maxcntD = max(max(den));
    if maxcntN > 2^16 || maxcntD > 2^16
        ovflw = 'yes';
        warning('Bin overflow warning');
    end
    
    %scaling the baselevel of the atom image
    numforscale=10;  %using a 10x10 pixel area to determine scaling fraction
    numleftup = sum(sum(num(1:numforscale,1:numforscale)));
    numleftdn = sum(sum(num(rows-numforscale:rows,1:numforscale)));
    numrightup = sum(sum(num(1:numforscale,cols-numforscale:cols)));
    numrightdn = sum(sum(num(rows-numforscale:rows,cols-numforscale:cols)));
    baselevelnum = numleftup + numleftdn + numrightup + numrightdn;
    
    denleftup = sum(sum(den(1:numforscale,1:numforscale)));
    denleftdn = sum(sum(den(rows-numforscale:rows,1:numforscale)));
    denrightup = sum(sum(den(1:numforscale,cols-numforscale:cols)));
    denrightdn = sum(sum(den(rows-numforscale:rows,cols-numforscale:cols)));
    baselevelden = denleftup + denleftdn + denrightup + denrightdn;
    blscale = (baselevelden/baselevelnum); %resulting baselevel scaling
    modnum=num*blscale;
    
    den(den==0) = 1; % Remove zeros from denominator
    
    ODraw = -log((num./den)');  % MATRIX ROTATED HERE
    ODrawmod = -log(modnum./den);
    ODraw(isinf(ODraw)) = 0;
    ODraw(isnan(ODraw)) = 0;
    ODraw = real(ODraw);
    
elseif type_img == 1 % Fluorescence
    ODraw = rawimages';         % MATRIX ROTATED HERE
    ODraw = ODraw - min(min(ODraw));
else
    error('Can''t determine imaging method');
end

% Filter image to find ROI
N = 30;  %amplitude of 2D gaussian filter
sig = 60;  %size of 2D gaussian filter
ODfilt = conv2(ODraw,gaussian2d(N,sig),'same'); % Low-pass filter image

yC = 1:expp.short_px_nb;   %HORIZONTAL DIMENSION OF EXPERIMENT
xC = 1:expp.long_px_nb;

autoROIok = false;

if autoROI % Find ROI and extract first approx. signal from filtered image
    % Determine pre-ROI for fitting (for initial guesses)
    mx = max(max(ODfilt));
    [I,J] = find(ODfilt == mx,1);
    wndw = 50;
    if ~(I - wndw < 1 || I + wndw > max(xC) || J - wndw < 1 || J + wndw > max(yC))
        autoROIok = true;
        
        ROI(1) = I(1);
        ROI(2) = J(1);
        
        % INTEGRATED ODfilt in diff. dimens. using pre-ROI as size (for initial guesses)
        yIntFilt = sum(ODfilt(:,ROI(2)-wndw:ROI(2)+wndw),2);
        xIntFilt = sum(ODfilt(ROI(1)-wndw:ROI(1)+wndw,:),1);
        
        % Least-squares fit of INTEGRATED ODfilt (initial guesses)
        ampGuessY = max(yIntFilt) - min(yIntFilt);
        ampGuessX = max(xIntFilt) - min(xIntFilt);
        options = optimset('Display','off','TolFun',1e-16,'TolX',1e-16);
        pYtmp = fminsearch(@(p)FitGauss(p,yC,yIntFilt),[ampGuessY ROI(1) 100 min(yIntFilt)],options);
        pXtmp = fminsearch(@(p)FitGauss(p,xC,xIntFilt'),[ampGuessX ROI(2) 100 min(xIntFilt)],options);
        
        % Determining autoROI based on size of ODfilt (for real fit)
        sigmas = 3;
        xWndw = round(sigmas*pXtmp(3));
        yWndw = round(sigmas*pYtmp(3));
        xROI = ROI(2)-xWndw:ROI(2)+xWndw;
        yROI = ROI(1)-yWndw:ROI(1)+yWndw;
        
        % Remove parts of the ROI if it goes outside of the photo.
        if min(xROI) < 1
            xROI = xROI(xROI > 0);
        end
        if min(yROI) < 1
            yROI = yROI(yROI > 0);
        end
        
        if max(yROI) > expp.short_px_nb
            yROI = yROI(yROI < expp.short_px_nb + 1);
        end
        if max(xROI) > expp.long_px_nb
            xROI = xROI(xROI < expp.long_px_nb + 1);
        end
        
        % INTEGRATED ODraw in diff. dimens. using autoROI as size (for real fit)
        yInt = sum(ODraw(yROI,xROI),2);
        xInt = sum(ODraw(yROI,xROI),1);
        
    else
        warndlg('Can''t auto-choose ROI: Estimated signal close to edge of image')
        autoROIok = false;
    end
end

if ~autoROIok  % Use user-input ROI parameters to get initial fitting parameters
    % User-defined ROI
    yROI = userROI(3):userROI(4);
    xROI = userROI(1):userROI(2);
    
    % INTEGRATED ODraw in diff. dimens. using autoROI as size (for real fit)
    yInt = sum(ODraw(yROI,xROI),2);
    xInt = sum(ODraw(yROI,xROI),1);
    
    % Amp. and sizes (initial guesses)
    pYtmp(1) = max(yInt) - min(yInt);
    pXtmp(1) = max(xInt) - min(xInt);
    pYtmp(2:3) = [mean(userROI(3:4)) 150];
    pXtmp(2:3) = [mean(userROI(1:2)) 150];
end

if doFitData % Real fit of ODraw using initial guesses
    options = optimset('Display','off','TolFun',1e-16,'TolX',1e-16);
    pY = fminsearch(@(p)FitGauss(p,yROI,yInt),[pYtmp(1:3) min(yInt)],options);
    pX = fminsearch(@(p)FitGauss(p,xROI,xInt'),[pXtmp(1:3) min(xInt)],options);
else
    pY = zeros(1,4);
    pX = pY;
end

if fit2d
    figure(100)
    subplot(2,3,1)
    surf(ODfilt,'Marker','.') %Plot filtered OD data
    cax=caxis;
    colorbar
    shading flat
    grid off
    title('Filtered ODraw')
    view(0,90);
    xlabel('vert')
    ylabel('horiz')
    
    % Binning ODraw for 2D fitting
    ODrawbin=BinImg2d(ODraw,bindim);  % NOTE - Rotated w/ respect to original image
    
    figure(100)
    subplot(2,3,4)
    surf(ODrawbin,'Marker','.') %Plot binned ODraw data
    cax=caxis;
    colorbar
    shading flat
    grid off
    title('Binned ODraw')
    view(0,90);
    xlabel('vert binned')
    ylabel('horiz binned')
    
    % Using ODrawbin for 2D fitting IC
    [horizbin,vertbin] = size(ODrawbin);  % y-dir. is HORIZ. exper. dir.; x-dir. is VERT
    amp2d=max(max(ODrawbin));
    [yoff2d,xoff2d] = find(ODrawbin == amp2d);
    sigx2d=vertbin/8;
    sigy2d=horizbin/8;
    guessave=3; %using a 3x3 binned pixel area to determine scaling fraction
    guessleftup = sum(sum(ODrawbin(1:guessave,1:guessave)))/guessave^2;
    guessleftdn = sum(sum(ODrawbin(horizbin-guessave:horizbin,1:guessave)))/guessave^2;
    guessrightup = sum(sum(ODrawbin(1:guessave,vertbin-guessave:vertbin)))/guessave^2;
    guessrightdn = sum(sum(ODrawbin(horizbin-guessave:horizbin,vertbin-guessave:vertbin)))/guessave^2;
    offset2d=(guessleftup + guessleftdn + guessrightup + guessrightdn)/4 ;
    slopex2d=((guessrightup + guessrightdn) - (guessleftup + guessleftdn))/(2*horizbin);
    slopey2d=((guessleftdn + guessrightdn) - (guessrightup + guessleftup))/(2*vertbin);
    IniCondNR=[offset2d,amp2d,xoff2d,yoff2d,sigx2d,sigy2d,slopex2d,slopey2d];
    
    % Creating matrices for 2D fitting with ODrawbin
    xvecbin=1:vertbin;
    yvecbin=1:horizbin;
    Nptsbin=vertbin*horizbin;
    [xdatabin,ydatabin]=meshgrid(xvecbin,yvecbin);  % xdata(ydata) have columns=length(xvec), rows=length(yvec)
    % xdata copies xvec along its rows; ydata copies yvec along its columns
    error2d(1:Nptsbin)=1;   % DUMMY
    npoints(1:Nptsbin)=1;   % DUMMY - Each entry is the number of points
    fitcoordmatrix=[xdatabin(:),ydatabin(:),error2d(:),npoints(:)];
    Z2dbin=ODrawbin(:);  %flattened(1D) vector of ODrawbin
    
    % Fitting - data is treated with equal weight
    opts=statset('TolFun',1e-6);
    [P2dbin,r2dbin,J2dbin]=nlinfit(fitcoordmatrix,Z2dbin(:),@gauss2derror,IniCondNR,opts);
    P2d = [P2dbin(1),P2dbin(2),P2dbin(3)*bindim(2),P2dbin(4)*bindim(1),P2dbin(5)*bindim(2),P2dbin(6)*bindim(1),P2dbin(7)/bindim(2),P2dbin(8)/bindim(1)];
    
    % Generate OD plot based on fit parameters
    [horiz,vert] = size(ODraw);  % y-dir. is HORIZ. exper. dir.; x-dir. is VERT
    xvec=1:vert;
    yvec=1:horiz;
    Npts=vert*horiz;
    [xdata,ydata]=meshgrid(xvec,yvec);
    ploterrorod(1:Npts)=1;
    plotnpoints(1:Npts)=1;
    plotcoordmatrix=[xdata(:),ydata(:),ploterrorod(:),plotnpoints(:)];
    
    ODcalcvec=gauss2derror(P2d,plotcoordmatrix);    %OD calculated from fit paramters
    ODcalc = reshape(ODcalcvec, size(xdata));
    
    ODcalcbinvec=gauss2derror(P2dbin,fitcoordmatrix);    %ODbin calculated from fit paramters
    ODcalcbin = reshape(ODcalcbinvec, size(xdatabin));
    
    ODcalcICvec=gauss2derror(IniCondNR,fitcoordmatrix);   %OD based on initial conditions
    ODcalcIC = reshape(ODcalcICvec, size(xdatabin));
    
    ODresiduals = ODfilt-ODcalc;
    
    subplot(2,3,5)
    surf(ODcalcIC,'Marker','.')  %Plot 2d fit using IC
    %caxis(cax)
    colorbar
    shading flat
    grid off
    view(0,90);
    xlabel('vert binned')
    ylabel('horiz binned')
    title('OD (binned) from initial conditions')
    hold on
    
    subplot(2,3,2)
    surf(ODcalc,'Marker','.')  %Plot 2d fit using fit param
    caxis(cax)
    colorbar
    shading flat
    grid off
    view(0,90);
    xlabel('vert')
    ylabel('horiz')
    title('OD from 2d fit')
    hold on
    
    subplot(2,3,6)
    surf(ODcalcbin,'Marker','.')  %Plot 2d fit using binned fit param
    caxis('auto')
    view(0,90);
    colorbar
    shading flat
    grid off
    title('OD (binned) from 2d fit param.')
    xlabel('vert binned')
    ylabel('horiz binned')
    hold off
    
    subplot(2,3,3)
    surf(ODresiduals,'Marker','.')  %Plot 2d fit residuals
    caxis('auto')
    view(0,90);
    colorbar
    shading flat
    grid off
    title('True residues (ODfilt-ODcalc)')
    xlabel('vert')
    ylabel('horiz')
    hold off
else
    P2d = [0,0,0,0,pX(3),pY(3),0,0];
end

% Give the width as the std deviation
if fit2d
    widthYm = (P2d(6))*expp.px_size*CCDbin/expp.mag;    %Y cloud size (m)
    widthXm = (P2d(5))*expp.px_size*CCDbin/expp.mag;    %X cloud size (m)
else
    widthYm = (pY(3))*expp.px_size*CCDbin/expp.mag;    %1D fits Y cloud size (m)
    widthXm = (pX(3))*expp.px_size*CCDbin/expp.mag;    %1D fits X cloud size (m)
end

if type_img == 1 % Fluorescence
    % Estimate number of atoms using camera
    if fit2d % if 2D fit has been done
        total_num_electrons = P2d(2)*2*pi*P2d(6)*P2d(5);
    else
        ampX = pX(1)/(sqrt(2*pi)*pY(3));
        ampY = pY(1)/(sqrt(2*pi)*pX(3));
        amp = mean([ampX ampY]);
        total_num_electrons = amp*2*pi*pX(3)*pY(3); % see p157 of Simon's lab book #1.
    end
    total_photons = total_num_electrons/(expp.coll_eff);
    numAtoms = 2*total_photons/(img_time*2*pi*cons.Gamma) % numAtoms = (num_photons/time) / (scattering rate)
elseif type_img == 3 % Absorption
    totalODraw = sum(sum(ODraw))*(expp.px_size*CCDbin/expp.mag)^2;
    numAtoms = totalODraw/AbsCross(0,0);
    if fit2d
        numAtomsFit = P2d(2)*2*pi*widthYm*widthXm/AbsCross(0,0)
    end
end

% Save desired information to handles structure to pass out of function:
handles.numAtoms = numAtoms;
handles.pX = pX;
handles.pY = pY;
if type_img == 3
    handles.totalOD = totalODraw;
    if fit2d
        handles.numAtomsFit = numAtomsFit;
    end
end

% Save desired information to HDF file:
SaveAnalysisHDF(hdfFilePath,'numAtoms',numAtoms);
if type_img == 3
    SaveAnalysisHDF(hdfFilePath,'totalOD',totalODraw);
    if fit2d
        SaveAnalysisHDF(hdfFilePath,'numAtomsFit',numAtomsFit);
    end
    SaveAnalysisHDF(hdfFilePath,'totalElec1',sum(sum(rawimages(:,:,1))));
    SaveAnalysisHDF(hdfFilePath,'totalElec2',sum(sum(rawimages(:,:,2))));
    SaveAnalysisHDF(hdfFilePath,'totalElec3',sum(sum(rawimages(:,:,3))));
end

SaveAnalysisHDF(hdfFilePath,'pX',pX);
SaveAnalysisHDF(hdfFilePath,'pY',pY);
SaveAnalysisHDF(hdfFilePath,'yData',yInt);
SaveAnalysisHDF(hdfFilePath,'xData',xInt);
if saveOD % Whether to save processed image or not (costly in terms of space)
   SaveAnalysisHDF(hdfFilePath,'OD',ODraw);
end


if calledFromGUI % If GUI called this function then update GUI display
    axes(handles.axes_processed)
    imagesc(ODraw)
    line([xROI(1) xROI(end)],[yROI(1) yROI(1)])
    line([xROI(1) xROI(1)],[yROI(1) yROI(end)])
    line([xROI(end) xROI(end)],[yROI(1) yROI(end)])
    line([xROI(1) xROI(end)],[yROI(end) yROI(end)])
    axis equal
    xlims = xlim;
    ylims = ylim;
    set(gca,'FontSize',10)
    
    axes(handles.axes_y_proj)
    plot(yInt,yROI,gauss(pY,yC),yC);
    set(gca,'YDir','reverse');
    ylim(ylims);
    set(gca,'FontSize',10)
    
    axes(handles.axes_x_proj)
    plot(xROI,xInt,xC,gauss(pX,xC));
    xlim(xlims);
    set(gca,'FontSize',10)
    
    if type_img == 3 % Abs. img.
        axes(handles.axes_raw_abs_atoms)
        imagesc(rawimages(:,:,1))
        axis off
        
        axes(handles.axes_raw_abs_no_atoms)
        imagesc(rawimages(:,:,2))
        axis off
        
        axes(handles.axes_raw_background)
        imagesc(rawimages(:,:,3))
        axis off
    else
        axes(handles.axes_raw_abs_atoms)
        imagesc(0)
        axis off
        
        axes(handles.axes_raw_abs_no_atoms)
        imagesc(0)
        axis off
        
        axes(handles.axes_raw_background)
        imagesc(0)
        axis off
    end
    
    set(handles.text_x_amplitude,'string',['X amp.: ' num2str(pX(1),'%.2e')]);
    set(handles.text_y_amplitude,'string',['Y amp.: ' num2str(pY(1),'%.2e')]);
    set(handles.text_x_sizepx,'string',['X width: ' num2str(pX(3),'%.2e')]);
    set(handles.text_y_sizepx,'string',['Y width: ' num2str(pY(3),'%.2e')]);
    set(handles.text_filename,'string',handles.dpath);
    set(handles.text_peak1d_position,'string',['peak(v,h) = (' num2str(round(pX(2))) ',' num2str(round(pY(2))) ')']);
    set(handles.text_roi,'string',['ROI(x,y) = (' num2str(xROI(1)) ',' num2str(xROI(end)) '),(' num2str(yROI(1)) ',' num2str(yROI(end)) ')']);
    set(handles.text_numAtoms,'string',['# Atoms: ' num2str(numAtoms,'%.2e')]);
    
    if type_img == 3
        set(handles.text_peak_pixel_count_signal,'string',['pk cnt(signal) = (' num2str(max(max(rawimages(:,:,1)))) ')']);
        set(handles.text_peak_pixel_count_bkgnd,'string',['pk cnt(bkgnd) = (' num2str(max(max(rawimages(:,:,2)))) ')']);
    end
elseif doDisplayImages
    figure(1)
    subplot(3,3,[1 2 4 5])
    imagesc(ODraw)
    line([xROI(1) xROI(end)],[yROI(1) yROI(1)])
    line([xROI(1) xROI(1)],[yROI(1) yROI(end)])
    line([xROI(end) xROI(end)],[yROI(1) yROI(end)])
    line([xROI(1) xROI(end)],[yROI(end) yROI(end)])
    axis equal
    xlims = xlim;
    ylims = ylim;
    set(gca,'FontSize',10)
    
    subplot(3,3,[3 6])
    plot(yInt,yROI,gauss(pY,yC),yC);
    set(gca,'YDir','reverse');
    ylim(ylims);
    set(gca,'FontSize',10)
    
    subplot(3,3,[7 8])
    plot(xROI,xInt,xC,gauss(pX,xC));
    xlim(xlims);
    set(gca,'FontSize',10)
end % end calledFromGUI